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Abstract 

The development of a thorough understanding of the mechanisms for vortex eruptions from viscous layers, 
which are believed to be associated with phenomena such as dynamic stall onset and transition, is crucial 
if accurate models of such phenomena are to be formulated. The development of such models may, in turn, 
allow for the possibility that such effects could be accounted for during the design of various aerodynamic 
devices such as wings, helicopter rotors and turbomachinery blading and thus lead to designs which are 
stall free or stall resistant and which have better stall-recovery properties. The present investigation is 
being conducted as part of an effort to develop analytical and numerical tools which can be used to help 
improve our understanding of the vortex-eruption mechanism at high Reynolds numbers. The addition of 
the normal-momentum equation to the classical unsteady boundary-layer equations is crucial, according 
to recent asymptotic analyses of the vortex-eruption problem, and is a key feature of the analyses being 
developed by the present authors. The purpose of this paper is, first, to describe departure solution behavior 
observed when using unsteady, streamline-curvature based solution procedures in which nontrivial transverse 
pressure gradient effects are included and, second, to show that special treatment of the time-derivative of 
the normal velocity is needed to eliminate the ill-posed solution behavior, which is observed when small 
spatial and temporal step sizes are used. 

Introduction 

A number of recent analytical studies have been directed towards understanding the fundamental physics 
associated with the development and subsequent eruption of concentrated regions of vorticity from the 
boundary layer (e.g., van Dommelen and Shen (1980), Elliott, et al (1983), Peridier, et al (1988) and Smith 
(1988)). This event is believed to be associated with well-known physical phenomena, such as the onset of 
airfoil dynamic stall and transition from laminar to turbulent flow. The cumulative observation of the above- 
mentioned and other studies seems to indicate that the classical boundary-layer equations are insufficient to 
completely describe the vortex eruption phenomenon. Even if strong viscous-inviscid interaction is allowed, 
it appears likely that normal pressure gradient effects must be accounted for in some form. The present 
paper describes part of an overall effort directed towards the development of unsteady analyses capable of 
addressing high Reynolds number flows in which normal pressure gradients are important and ultimately, it 
is hoped, where vortex eruptions occur, as well. 

The important work of van Dommelen and Shen (1980) first documented the existence of a finite-time 
singularity in the solution of the non-interacting, classical, unsteady boundary-layer equations for flow past 
an impulsively started circular cylinder. Later work (e.g., Peridier, et al (1988)) showed that, if the boundary- 
layer equations are allowed to interact with the inviscid flow, the van Dommelen and Shen singularity can 
be bypassed. However, another finite-time singularity, which cannot be removed through interaction alone, 
arises shortly thereafter. The recent asymptotic study of Smith (1988) indicates that normal pressure gradient 
effects must be included in order to avoid the latter (interactive) finite-time singularity. This provides the 
motivation for the present work, which addresses the issue of how various terms in the unsteady normal- 
momentum equation must be numerically treated within the framework of a globally it erated space- and time- 
marching solution algorithm, in order to avoid ill-posed behavior of the solutions. Globally iterated solution 
algorithms for steady and unsteady flows have been employed extensively in both asymptotic (i.e., infinite 
Reynolds number) and finite Reynolds number investigations, where they have proven to be computationally 
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efficient. Streamline curvature techniques have been demonstrated in a number of studies, for example, see 
Smith ci al (1984), Presz (1985), Rothmayer (1989, 1990) and Power (1990) - many other examples can be 
found in the literature. Thus, this approach is being pursued here in the hope that the use of a boundary- 
layer like numerical approach will lead to a computationally efficient technique to solve the vortex-eruption 
problem. 


External Flow Analysis - Flow Past a Flat Plate 

The particular concern of this paper is an issue which arose while the first author was developing a 
numerical solution scheme for the equation set consisting of the classical, unsteady, incompressible boundary- 
layer equations supplemented with the inviscid form of the normal-momentum equation. Note that these 
equations are essentially identical to the leading-order terms in the Incompressible form of the “Thin-Layer” 
^.avier-Stokes equations, and will therefore be referred to here as the ITLNS equations, for convenience. For 
two-dimensional flows these equations are given by 
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where u and v are the velocity components in the x- and y- direct ions, respectively, with x oriented tangent 
to the body surface and y normal to it, and P is the static pressure. In this section, the body is assumed to 
be a semi-infinite flat plate, so that x and y are Cartesian coordinates. Standard low-speed, external flow 
nondimensionalizations have been used and the y-coordinate and v-velocity component have been scaled with 
the square root of the Reynolds number as follows: u = u* /U^, v = v*y/Re/U^ 0 , P = (P* — Poo )/p* Uoo j 
x — x* /L* e j, y = y’VRe/ L‘ e j and t = t* /{L* e j /U^). Asterisks denote dimensional quantities, the subscript 
oo denotes a quantity evaluated in the uniform far field upstream flow and Re is the Reynolds number defined 
as — UZaL’ej/v*, where v* is the upstream value of the kinematic viscosity which, along with the density 
p* , is assumed to be constant. 

The ITLNS equations can be solved in the primitive variable form given above, or they can be solved after 
transforming to Gortler variables — the latter approach has been used in this study. However, to simplify the 
present discussion, the primitive variable form of these equations will be considered, after using the stream 
function V’> defined by the relations u = dip/dy and v = —dipjdx, to replace v. Substituting for v in Eq. (3) 
yields 

d 2 ip d 2 ip dip d 2 ip _ n ®P 

dxdt U dx 2 dx dxdy dy 


(4) 


The boundary conditions on the surface (y = 0) are the no-slip, zero injection conditions. u(x,0) 0 

and 4<(x, 0) = 0 for t > 0. At the outer edge of the boundary layer (y -*• oo) the edge condition on u is 
lim u(x,y;t) — * U e (x,t) and the pressure P satisfies the unsteady Bernoulli relation. The equations can 

either be solved in direct mode (edge conditions specified) or inverse mode (displacement thickness specified). 
The additional boundary condition needed for the latter is given by the following relationship between the 
edge value of ip and the displacement thickness 6 * , 


4’(x,y e ;i) = U e (x,t)(y e -6*(x,t)) 


(5) 


where y e is the value of y at the boundary-layer edge and U t is treated as an unknown. In addition, both 
upstream and downstream boundary conditions are needed - the latter is required because the introduction 
of the inviscid normal-momentum equation makes the in vise id form of the governing equations equivalent 
to the unsteady incompressible Euler equations, which are elliptic-like in space at any time f, as explicitly 
evidenced by the presence of the ip IX term in Eq. (4), which represents the streamline curvature. For the 
simple problem to be considered here, the upstream and downstream profiles for u and V' are assumed to 
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correspond to the Blasius (flat-plate) profile. Finally, the initial condition which is needed at time t — 0 is 
assumed to correspond to the Blasius solution along the entire plate surface. 

The ITLNS equations are numerically solved using fully implicit discretizations which are first-order 
accurate for the x- and ^-derivatives and second-order accurate for the y- derivatives. All y - derivatives are 
central differenced and all x- and ^-derivatives are backward differenced with respect to the solution point 
(at (x,-, yj , t n )) with the exception of the 4 y xx term in the normal-momentum equation, which is centra] 
differenced, thus introducing unknown information from the upstream point at x J+I . This information 
is obtained by initially guessing ^i+i> using the value from the previous time step, and then performing 
multiple, global spatial sweeps at each time level until the ^-field converges. An acceleration scheme has 
been used to improve the convergence properties of this procedure but will not be discussed here, however, 
as it is not relevant to the present focus. Initially the authors believed that central differencing of the $ xx 
term would properly (and fully) account for the elliptic-like nature of the governing equations at each time 
step, as demonstrated in applications of similar approaches to steady flows (e.g., see Presz (1985)). However, 
this was not found to be the case, as will be discussed below. 

The issue which is of concern here arises when attempting to solve the above system of equations for 
a specified time-dependent displacement-thickness distribution (i.e., the inverse method). Consider the 
simplest possible case, where the displacement thickness is assumed to be that for a steady flat-plate flow for 
all time t > 0. Thus, the inverse solution procedure should yield the Blasius solution at all values of x and t , 
with a small perturbation (depending on the specified value of the Reynolds number) due to normal pressure 
gradient effects. This has been found to be the case here when the spatial and temporal step sizes Ax and A t } 
respectively, are not chosen to be too small. However, as Ax and A t are decreased, it has been observed 
that, during the first time step, the solution departs from its anticipated behavior in a manner reminiscint 
of that observed when attempting to solve a boundary-value problem using an initial-value technique. This 
occurs despite the use of central differencing for the term. That is, for a fixed, constant spatial stepsize 
Ax below some minimum value, there appears to be a minimum temporal step size At , below which the 
space-marching solution behaves as if it is ill-posed with respect to x. 

Examples of the departure solutions are shown in Figs. 1A and IB, where the skin-friction coefficient, 
Q = 2r* V^e/ 'p*U^ , and wall pressure P w , respectively, are plotted as functions of distance along the plate 
for a Reynolds number of 1 x 10 6 based on a reference length L ref - 1. This case was calculated starting 
at x = 1.0 with a fixed value of Ax = 0.001 and three different values of At, namely 0.0010, 0.0009 and 
0.0001. This value of Ax is below the minimum for which departures have been observed, and the three 
values of At are near the boundary between departure-free solutions and departing solutions. Note that the 
solution goes from being well-behaved at the largest value of At to growing in an oscillatory exponential 
manner for the middle value to monotonic exponential growth for the smallest value. Similar departure 
behavior was subsequently observed in the solutions obtained from a different numerical code which uses a 
similar streamline-curvature based technique to solve the full Navier-Stokes equations, for interna] flows, as 
discussed in the next section. 

Before continuing, it should be noted that a consistency check on the finite- difference form of Eq. (4) was 
carried out. The equations were found to be consistent in the sense that as Ax and A? are independently 
reduced to zero, Eq. (4) is identically recovered. Thus, the possibility that truncation errors associated with 
the discretized equations have changed the mathematical character of the governing differential equations 
has been eliminated as a possible source for the branching behavior. 

Calculations to establish the departure boundary” for three different Reynolds numbers were performed, 
namely Re — lx 10 6 , 1 x 10 7 and 1 x 10 8 . The results are consistent, and indicate that the value of At 
for which the solution crosses over from “departure-free” to “departing” is a function of Reynolds number, 
as might be anticipated. The fact that there is a minimum Ax above which solutions remain departure-free 
for any value of At is not surprising - once Ax becomes large enough, the solution probably oversteps the 
streamwise length scale of the physical mechanism governing the behavior. 

The branching behavior of the small step size solutions of the ITLNS equations obtained using the present 
numerical solution procedure has been examined in detail by the authors in an effort to understand its source. 
As a result of this investigation, the term responsible for the departure solutions has been found to be the u t 
term appearing in the normal-momentum equation. That is, it has been found that, regardless of the values 
used for Ax and At , if the Vt term is neglected , then the solution will never exhibit the departure behavior 
described above, so long as the ip xx term is not backward differenced, but is instead central differenced. 
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Further, this behavior is found to be independent of the numerical treatment used for the convective terms 
in the normal-momentum equation. 

Recall that in the present study, the v t term has been written using the stream function definition, i.e., 
v t = -ip xt . In the numerical algorithm described above, this term was discretized using a backward difference 
for the a*- derivative, leading to the following form at all y-locations, where the subscript i and superscript 
n denote the x-index and the f-index, respectively, with (/, n) at the current station in both space and time 
and Ax assumed to be uniform: 


4>xt 


At Ax 


[W-eO-W'-C"! 1 )] 


( 6 ) 


The apparent ill-posedness exhibited by the numerical solutions, and the fact that the v t term has been 
found to be responsible, suggests an alternative discretization for the V’xt term wherein a forward difference 
in x with respect to the solution station is used. That is, tl> £t is discretized in the form 




I 
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With this modification, the ill-posed behavior which is observed when a backward difference is used no 
longer arises, regardless of the values of Ax and At that are specified. A case for which violent branching 
occurs when using a backward difference for the x-derivative in Vvt, which was presented in Fig. 1, has been 
recalculated using a forward difference with the same spatial and temporal stepsizes, i.e., Ax = 0.001 and 
At = 0.0001. The resultant solution is departure free and virtually identical to the backward-difference 
solution for At = 0.001, which did not branch because the value of Ax was too large. 

The precise reason that forward-differencing of the x-derivative in the ip st term is needed to prevent 
the ill-posed behavior of the small step-size solutions of these equations is unknown at the present time. 
One possiblity is that the responsible mechanism is somehow related to physical boundary-layer instability 
mechanisms. Another possibility is that the mechanism is purely inviscid in nature, as is the case leading 
to the requirement for central differencing of the il ] X x term. Both of these possiblities are currently under 
investigation, and it is hoped that this issue will be resolved in the near future. 


Internal Flow Analysis - Pulsatile Flow Through a Channel 

Here the governing equations are taken to be the unsteady Navier-Stokes equations. This set of equations 
is non-dimensionalized by the method of Smith (1976). The final equations are: conservation of mass, 

+ v y = 0 ; (8) 


conservation of x-momentum, 

U t + Re(uu x + VUy) = —Par + Wjx + Uyy '■> 


(9) 


and conservation of y-momentum, 

V t + Re[u(v x — Vt) + ttVy] = —Py + + v yy • (^) 

The v T term is a pseudo-time derivative introduced to accelerate convergence of^the global-iteration scheme 
at each time-level t. The Reynolds number, Re, is defined here by Re = L* g* /p*v* , where L * is the 
dimensional channel width, p* is the density, u* is the kinematic viscosity and -g w is the local applied 
pressure gradient driving the basic flow. This set of equations is solved in a two-dimensional channel where 
the upstream flow is a pulsatile Poiseuille flow driven by the pressure gradient 

a p 

— (-oo,j/) = -1 +7o cos j3l . (11) 

OX 

The corresponding velocity profile is 

u(-oo.y) = u o + u o > 0^) 


where 


K = ly(l-y) 


iS = Cie v^-y + C 2 e ^y_J | 


(13) 
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and v(-oo,y) = 0. No-slip boundary conditions are applied along the upper and lower walls and the 
downstream boundary conditions are n x = 0 and v x = 0. 

For the case where the streamwise length scale is the same order as the channel width, the minimal system 
of equations needed to reproduce the unsteady asymptotic structure of Smith (1976) is the parabolized 
Navier-Stokes (PNS) equations (Eqs. (8)— ( 10) with all underlined terms neglected). Bearing this in mind, 
a Navier-Stokes algorithm is formulated by first developing a PNS algorithm and then iterating on the 
additional (underlined) terms. It is well-known that, when solving the steady PNS equations, departure 
solutions can be suppressed either by neglecting the streamwise pressure-gradient term P X) or the streamline- 
curvature term (the uv x term, see Rothmayer (1989, 1990)). In this study, the streamline-curvature term 
uv x is used to suppress departure solutions in the steady single-pass version of the algorithm, and also to 
provide a mechanism for upstream propagation of information in the steady and unsteady global-iteration 
algorithms. 

The streamline-curvature term is treated as a known source term and is forward-differenced in space 
relative to the current solution point. The pseudo-unsteady term uvx (introduced to accelerate convergence 
as in Davis (1984) and Barnett and Davis (1986)) is included implicitly using a standard backward difference 
in time (see Rothmayer (1989, 1990)). Note that two unsteady effects are present - the real unsteady terms 
u t and v t and the pseudo-unsteady term Vj , The latter is driven to zero during multiple sweeps through the 
solution domain at each real time level t. 

The algorithm described above, without the addition of the Navier-Stokes terms (i.e., underlined terms 
in Eqs. (8)-(10)) is similar to that described in the section on external flows, where the ITLNS equations 
are solved. As with that algorithm, a number of implicit/explicit PNS-like algorithms were tested. The 
full Navier-Stokes version of the internal-flow solution technique treats the underlined terms as source terms 
calculated from the solution at the previous iteration, although algorithms with implicit treatment of v t were 
also tested. The reader is referred to Rothmayer (1990) for further details. 

The above-described Navier-Stokes solution algorithm has been used to solve for the flow through a 
flat channel with the pulsatile pressure gradient given by Eq. (11) and the upstream boundary conditions 
given by Eqs. (12-14). Fig. 2A shows a comparison between the wall shear stress computed using the 
present analysis at a downstream location along with that given by the analytical Poiseuille solution - the 
agreement is excellent. The departure solutions, to be discussed next, were triggered by introducing a very 
small indentation in the channel wall (typical height h - lx 10~ 6 ). 

As with the previously described external-flow analysis, it has been found that the present solution 
algorithm experiences departure-solution behavior when a minimum spatial/temporal step-size restriction is 
violated, with solutions like that shown in Fig. 1. Figure 2B shows how the minimum allowable time step, 
A/j, changes with varying streamwise step size, A#, for a Reynolds number of 10 million. For a streamwise 
step size above a critical value ( Ax % 0.207) the numerical scheme is free of departure solutions for all values 
of At examined. As found with the external-flow analysis, neglecting the v t term leads to departure-free 
solutions for all values of Ax and At. 

A similar departure-solution behavior was also observed when solving the PNS equations numerically. As 
with the Navier-Stokes algorithm, the unsteady PNS method displayed the spatial/temporal step-size con- 
straints. These departure solutions could again be eliminated by neglecting the v t term. These observations 
hold even if the normal-momentum equation is reduced to the very simple form v t = -P y . 

In the external-flow analysis, branching of the numerical solution was suppressed by spatially forward- 
differencing the v t term, after re-expressing in terms of the stream function. A similar approach was at- 
tempted in the internal-flow analysis. The —4 } xt term was forward and backward differenced in space, and 
treated both implicitly and explicitly in both cases, in an attempt to eliminate the departure solutions. Of 
the four methods, the backward-differenced explicit method had the least severe time-step restriction for a 
given value of Ax. However, for sufficiently small time steps, all four algorithms exhibited the departure- 
solution behavior previously discussed. This is in contrast with the external-flow analysis, where branching 
could be completely eliminated by spatial forward-differencing of the term. It is not clear at this time 
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why this difference between the internal- and external-flow analyses exists. Fortunately, the pseudo time 
step AT can h>e optimized so that the time scales for the observed departure solutions fall T>elow the scales of 
the Navier-Stokes regime (which constitute a likely absolute lower bound on the temporal time step needed 
for practical calculations). However, if a judicious choice of AT is not made, then departure Solutions may 
be encountered even at the large values of At associated with the interactive boundary-layer regime of Smith 
(1976) (see Fig. 2C). 

Anothet interesting, and currently unexplained, phenomena observed in both the internal- and external- 
flow analyses is that departure solutions can be eliminated by neglecting the v* Term in certain local transverse 
regions, while retaining this term outside of those regions. The five points labeled Al, Xi andl3l through 
B3 on Fig. 2B are points at which the solution has been stabilized by neglecting Vt in the various regions 
indicated. Fig. 2D shows the location of these regions for each point. It can be seen from This group of 
figures that the size and location of these regions are sensitively dependent on both the spatial and temporal 
stepsizes. For all values of A.r there appears to exist a At below which v t must be neglected across the entire 
channel to ensure departure-free solutions. 


Concluding Remarks 


The objective of this paper has been, first, to describe departure solution behavior observed when using 
streamline-curvature based solution procedures which are being developed to study high ‘Reynolds number 
vortex-eruption phenomena, second, to indicate the responsible term in the governing equations and, finally, 
to show how the departure solutions can be eliminated. We have shown that the time-derivative of the 
normal velocity, v t = —tyxu appearing in the normal-moment um equations, is responsible for the branching 
behavior, which only occurs for small spatial and temporal step sizes. The ill-posed behavior has been 
eliminated in the external-flow analysis by forward differencing the spatial-derivative appearing in the — ij> rt 
term. It should be noted that the step sizes for which the observed ill-posed behavior arises turn out to 
be within the range needed to capture many important unsteady viscotis-inviscid interaction phenomena, 
such as dynamic stall onset. Therefore, this mechanism should not be ignored if accuraTe solutions are to be 
obtained. The implication is that special differencing procedures may be needed to properly account for the 
mechanisms responsible for the elliptic-like character of the governing equations at each time level of a time- 
marching algorithm, possibly even for non-streamiine-curvature techniques. A more complete description 
of the responsible mechanism is currently being pursued by the authors, and will be reported when it is 
available. 
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Distance along plate, x 

Fig. 1A. Wall pressure along plate for three values of At, Ax=0.001. 



Fig. IB. Skin-friction coefficient along plate for three values of At, Ax=0.001 
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Fig. 2A. Comparison of analytical and calculated wall shear with the solution driven 
pulsatile pressure gradient. 
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Fig. 2C. Variation of minimal time step constraint with changing pseudo time step 
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Fig. 2D. Illustration of regions where neglecting V t suppresses branching. 
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